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ABSTRACT 


We address problem of autonomous cryogenic management 
of loading operations on the ground and in space. As a step 
towards solution of this problem we develop a probabilistic 
framework for inferring correlations parameters of two-fluid 
cryogenic flow. The simulation of two-phase cryogenic flow 
is performed using nearly-implicit scheme. A concise set of 
cryogenic correlations is introduced. The proposed approach 
is applied to an analysis of the cryogenic flow in experimen- 
tal Propellant Loading System built at NASA KSC. An effi- 
cient simultaneous optimization of a large number of model 
parameters is demonstrated and a good agreement with the 
experimental data is obtained. 


1. INTRODUCTION 


NASA program of advanced cryogenic loading operations is 
focused on development of technologies enabling effective 
space exploration (Chato, 2008; Notardonato, 2012; Konishi 
& Mudawar, 2015a). These technologies require models that 
can recognize and predict cryogenic fluid regimes without hu- 
man interaction. The accuracy of these predictions depend 
largely on the functional form and parameterization of cor- 
relations for the flow patterns, pressure losses, and the heat 
transfer in cryogenic fluids. However, knowledge of these 
correlations remains limited and their learning present sig- 
nificant challenge to scientists and engineers (Tong & Tang, 
1997; Faghri & Zhang, 2006; Ghiaasiaan, 2007; Konishi & 
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Mudawar, 2015b) This problem becomes even more chal- 
lenging for applications in micro-gravity (S. R. Darr et al., 
2016a, 2016b). 


Research conducted to address this and other critical issues of 
space exploration resulted in a number of important achieve- 
ments, see e.g. (Konishi & Mudawar, 2015a; S. R. Darr et 
al., 2015; Hartwig & Vera, 2015; Hedayatpour, Antar, & 
Kawaji, 1990; Kawaji, 1996; A. Majumdar & Steadman, 
2003; A. K. Majumdar & Ravindran, 2010; Yuan, Ji, & 
Chung, 2007, 2009; Chung & Yuan, 2015; S. Darr et al., 
2016; S. R. Darr et al., 2016a, 2016b) and references therein. 
Specifically, a number of optimization techniques have be- 
come available for analysis of the model parameters and data 
correlations (Cullimore, 1998). 


However, small time steps and instabilities of numerical 
schemes impose substantial limitations on the speed of the 
two-fluid cryogenic codes and possibility of their on-line ap- 
plications in autonomous systems. As a result accurate pre- 
dictions of transient cryogenic flows and optimization of the 
parameters of cryogenic systems remain a challenge (Hartwig 
& Vera, 2015; Cullimore, 1998; Konishi & Mudawar, 2015a; 
S.R. Darr et al., 2015). 


Here we report on progress in development of probabilis- 
tic framework capable of efficient multi-parametric inference 
of cryogenic correlations and accurate predictions of experi- 
mental time-series data. 


The paper is organized as follows. In the following section we 
briefly describe model equations and method of their integra- 
tion. In Sec. 3 we discuss correlation relations for two-phase 
cryogenic fluid. In Sec. 4 we introduce probabilistic frame- 
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work for multi-parametric inference of the model parameters. 
An application of the inferential framework to the analysis of 
cryogenic flow in experimental Propellant Loading System is 
discussed in Sec. 5. Finally, in Conclusions we summarize 
the obtained results and discuss future applications. 


2. MODEL 
2.1. Model equations 


The system is modeled using six conservation equa- 
tions for gas/vapor and liquid phases in a set of con- 
trol volumes (D. G. Luchinsky, Smelyanskiy, & Brown, 
2014a, 2014c, 2014b; D. G. Luchinsky et al., 2016), see 
also (Luchinskiy, Ponizovskaya-Devine, Hafiychuk, et al., 
2015; D. G. Luchinsky et al., 2015). The set of equations can 
be written as following 


9 | cow | + 2 ope? 
ape |, apeu 
0 
2 0 (1) 
p (sa, — (au) ;) 
sj; T 
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Here 7 takes values g for the gas and / for the liquid and dy- 
namical variables in this system are: void fraction a (volume 
conservation requires ag + a; = 1), pressure p, phasic den- 
sities p,(y), velocities ug), and thermal energies e,(1), 8; is 
negative sign for the liquid and positive for the gas. The sys- 
tem (1) is closed by constitutive relations and equations of 
states implemented in the form of tables. The values of heat 
(q) and mass (TL) fluxes at the wall and interface (subscripts w 
and 7) are given by correlation relations. 


The system (1) is coupled to the energy conservation equation 
for the pipe walls 


Pwlwdy ee = Pigg (Ty Tw) 


2 
t+hwi (Ty = Tw) + hamb (Liane = Tis) : ( ) 


Here T,, is the pipe wall temperature, py, Cy, and d, are 
density, specific heat, and thickness of the pipe wall, h is the 
heat transfer coefficient corresponding to the ambient (hamp) 
and internal heat flowing to the wall from the gas (h,,,) and 
liquid (h,1) phases. 


The total mass transfer I‘, in equations (1) is the sum of the 
mass transfer at the wall (w) and at the interface (2) 


Pg =Dug + Pig, 


where 
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The heat fluxes at the wall and at the interface are defined as 
follows 


dwg = hwg (Tw — Tg) : 
dui =o hwi (Tw a Ti); 


Here T/,, is the saturation temperature of the liquid. 


dig = hig (Ths -_ Ty) ; 
Gat = hu (Tis — Th) - 


2.2. Algorithm 


The integration of the model (1) - (2) was performed in two 
steps following closely the nearly-implicit method described 
in RELAP5-3D (RELAP5-3D-I, 2012), see (D. G. Luchinsky 
et al., 2014c; Luchinskiy, Ponizovskaya-Devine, Hafiychuk, 
et al., 2015) for the details. The first step of the algorithm can 
be briefly summarized as follows: (i) Solve expanded equa- 
tion with respect to pressure in terms of new velocities; (ii) 
Solve momenta equations written in the form of block tri- 
diagonal matrix for the new velocities; (iii) Find new pres- 
sure; (iv) Find provisional values for energies and void frac- 
tions using expanded equations; (v) Find provisional values 
of mass fluxes and heat transfer coefficients using provisional 
values of temperatures obtained. 


At the second step new values of the densities, void fractions, 
and energies are found by solving the unexpanded conserva- 
tion equations for the phasic masses and energies using pro- 
visional values for the heat and mass fluxes in source terms. 
The solution is reduced to independent solution of four tri- 
diagonal matrices. The values of pressure and velocities in 
these matrices are taken at the new time step. 


The resulting computational scheme is efficient and fast and 
can integrate 1000 sec of real time chilldown in a few seconds 
of computational time. For a model consisting of N control 
volumes it involves inversion of N 4 x 4 matrices, solution of 
2 x (N — 1) tree-block-diagonal matrix equation, solution of 
four N x N tridiagonal matrix equations, and N x m explicit 
computations. 


The mitigation of known stability issues of this algorithm is 
discussed in (D. G. Luchinsky et al., 2014a, 2014c, 2014b, 
2016; Luchinskiy, Ponizovskaya-Devine, Hafiychuk, et al., 
2015; D. G. Luchinsky et al., 2015) 


3. CONSTITUTIVE RELATIONS 


The model equations introduced above have to be closed 
using the equations of state and the constitutive relations. 
The equations of state can be included into the model in the 
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form of NIST tables (Lemmon, Huber, & McLinden, 2013; 
D. G. Luchinsky et al., 2014c). Functional and parametric 
form of constitutive relations, on the other hand, represent 
one of the main source of uncertainties in the model. In this 
work we use constitutive relations (D. G. Luchinsky et al., 
2016) based on the recognition of the flow boiling regimes 
and for each regime we define interphase friction, the coef- 
ficient of virtual mass, wall friction, wall heat transfer, inter- 
phase heat and mass transfer. 


At each time step the boundaries between flow boiling 
regimes are determined first, then the flow correlations are 
defined for each flow regime. 


The boundaries between flow regimes are estimated using 
simplified Wojtan et al (Wojtan, Ursenbacher, & Thome, 
2005) map. The map was simplified in two ways. First we 
used only three transition boundaries. Next, we estimated the 
location of these boundaries in the coordinates of mass flow 
rate (7m) and vapor quality (v) using original expressions. Fi- 
nally, we approximated the location of these boundaries us- 
ing low-dimensional polynomials and used polynomial coef- 
ficients as fitting parameters. 


The rationale behind approximate description of the bound- 
aries of flow regimes is twofold. It is known (Jackson, 2006) 
that Wojtan maps can only be considered as an approximation 
to the flow regimes for cryogenic fluids. In the experiments 
performed at KSC the flow regimes were not established ex- 
perimentally and therefore cannot be validated. 


The following transitions between flow regimes were in- 
cluded into correlations: (i) Stratified- Wavy-to-stratified tran- 
sition; (i) Stratified-Wavy-to-annular-intermittent transition; 
and (iii) Dryout transition. 


3.1. Heat and mass transfer 


In the current work we are interested in relatively low mass 
fluxes, G < 600kg/m?/s. In this limit correlations for the 
heat flux are often based on the multiplicative or additive cor- 
rections to the values obtained for pool boiling (Franchello, 
1993; Seader, Miller, , & Kalvinskas, 1965; Griffith, 1957; 
Zuber & Tribus, 1958; Kutateladze, 1959). The following 
heat transfer mechanisms are included in the analysis: (i) con- 
vection, (ii) nucleate boiling, (ii) transition boiling, (iv) film 
boiling, and (v) transition to dryout regime. 


3.1.1. Convective Heat Transfer 


Convective heat transfer in horizontal pipes distinguishes four 
flow regimes: forced laminar, forced turbulent, natural lami- 
nar, and natural turbulent convection. The corresponding cor- 
relations for the convective heat transfer can be taken in the 
form e.g. (TRACE-V5.0, 2007; RELAP5-3D-IV, 2012; Nel- 


lis & Klein, 2009; Holman, 1989) 


4.36, Forced-Laminar; 
K 0.023 - Re®-8 Pr®-+, Forced-Turbulent; 
D;,, \ 0.1-(Gr- Pr)*/, Natural-Laminar; 
0.59 - (Gr - Pr)'/4, Natural-Turbulent. 


heb = (3) 


2 3 
_— He. _ p98 (Tw=Tig))D 
Here Pr = “” and Gr = — 


= 7 are Prandtl 
and Grashof numbers respectively, Gr is the coefficient of 
thermal expansion, and Dy» is the hydraulic diameter. To 
guarantee a smooth transition between various regimes the 
maximum value of h,, is taken as the value for the convec- 


tive heat transfer. 


We note that convective heat transfer in the stratified flow 
does not significantly affect the chilldown process, because 
the fluid temperature in this regime is close to (or lower than) 
saturation temperature T’,. The first critical temperature that 
defines the shape of the boiling curve and influences the chill- 
down corresponds to the onset of nucleation boiling To,,»5. 


3.1.2. Onset of nucleate boiling 


The correlations for temperature T;,,,, and the correspond- 
ing heat flux gon» for onset of nucleate boiling can be writ- 
ten (Sato & Matsumura, 1964; Frost & Dzakowic, 1967; Ghi- 
aasiaan, 2007; Huang, 2009) as 


IAT. 
Ton = nae(14 es ‘). (4) 
: B 2 
donb = Pr2 AT: 4 = heo(Tonb _ T;) (5) 


2 
where B = Sate"! P= ASEM AT. 4, = Ton — Ts is the 
wall superheat, and AT,,,, = T, — Tj; is liquid subcooling 
temperature. The convective heat transfer coefficient is given 


by (3). 


When the wall temperature is increased the heat flux to the 
wall may increase by more than an order of magnitude until 
the heat flux approaches its critical value qn f. 


3.1.3. Critical heat flux 


The values of critical heat flux q.,f and the correspond- 
ing critical wall superheat T(.,¢ are crucial for predicting 
chilldown and dryout phenomena in non-equilibriums flows. 
In nuclear reactor codes. (RELAP5-3D-IV, 2012; TRACE- 
V5.0, 2007) these values are determined using look-up tables 
based on extensive experimental measurements obtained un- 
der various flow conditions. For cryogenic fluids experimen- 
tal data remain sparse and values of q.,7 and Tiny are of- 
ten estimated using mechanistic models (Tong & Tang, 1997; 
Ghiaasiaan, 2007; Seader et al., 1965; Crowley & Izenson, 
1989; Konishi & Mudawar, 2015b). 
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The temperature T(,), ¢ for the critical heat flux was estimated 
in this work using approach proposed by Theler (Theler & 
Freis, 2011) 


Ts 


; (6) 
1 — 42 }0g(2k, + 1) 


Tenf = 


where kz, is the isoentropic expansion factor that for ideal di- 
atomic gases is 7/2 and R, is the specific gas constant. 


In boiling flows critical heat flux corrections have to take into 
account the dependence of the heat flux on the void frac- 
tion, velocity, and sub-cooling of the flow. In this work we 
used correlations proposed by Griffith et all for cryogenic 
flows (Franchello, 1993; Seader et al., 1965) 


. CAT sx 
Qchf = ch f,0(Qer _ a) (1 + ay (sent) (7) 
Pghig 
Rempel vt 
ey eee (Aue | 
Pahig 


where a, is the critical value of the void fraction and a; 
are constants, e.g. a; = 0.0144, ag = 10-8, a3 = 0.5 x 
1073 (Griffith, 1957), and a,, = 0.96 (Franchello, 1993) for 
water. 


The term dcp f,o in the eq. (7) corresponds to the pool boiling 
correlations, which were estimated in the present work using 
Zuber (Zuber & Tribus, 1958) correlations 
1/2 
) . (8) 


og(pi— fo) ) io ( pl 
PG pi+ Pg 
In numerical tests, we often used a simplified expression for 
eq. (7), cf (Franchello, 1993; Iloeje, Plummer, Rohsenow, & 
Griffith, 1982) 


: T 
dchf,0 = 54 ltlaPo ( 


Ach f = Ich f,0 “ay° (Qer = a)? (1 + a3G®*), (9) 


where typical values of parameters used in simulations are 
a1=1.0, A,=0.96, a2=2.0, ag=0.16, and a4=0.2. 


When wall superheat exceeds AT.,¢ = Teng — Ts, the tran- 
sition boiling begins and the heat flux to the wall decreases 
sharply as a function of the wall temperature until the latter 
reaches minimum film boiling temperature T,,, fp. 


3.1.4. Minimum film boiling 


In the film boiling regime the fluid flow is completely sepa- 
rated from the wall by the vapor film. The minimum value 
of the wall superheat AT; ¢o = Imp — Ts corresponding 
to this regime was estimated by Berenson as (Carbajo, 1985; 


Berenson, 1961) 


AT mfo,0 = 0.127 Pete x 
g(pi—Pg) are o . Hg ue (0) 
| pita lata | lata 


Tloeje (Franchello, 1993; Iloeje et al., 1982) has corrected 
Berenson equation to take into account the dependence of the 
AT fp on the quality and mass flux of the boiling flows in 
the form 


AT m fo — a ATn Fb,o(1 _ coX 63) (1 + caG®), (11) 


where X, is the equilibrium quality, G is liquid mass flux 
and a, are constants, e.g. a; = 0.0144, ag = 10-©, ag = 
0.5 x 1073 (Griffith, 1957), and a, = 0.96 (Franchello, 
1993) for water. 


The heat flux in the film boiling flow can be chosen 
following e.g. recommendations of Groeneveld and 
Rousseau (Groeneveld & Rousseau, 1983). In this work 
the heat flux to the wall in the film boiling regime was taken 
in the form of Bromley correlations 

IPgks (pt — Pg) higepg ca 


hor =C- ; 
: Dil — To Pry 


(12) 


corrected using Iloeje-type correlations (Franchello, 1993; 
Tloeje et al., 1982) 


hyo = chp, (1 -— @XS)(14+ 4G) (13) 


Typical values of the parameters used in simulations are the 
following: (i) c, = 2.0; (ii) cg = 1.04; (iii) cz = 2.0; Gv) 
c4 = 0.2; (v) cs = 0.1. 


The minimum film boiling heat flux can now be defined as 


dmb = hypo ATmfo- (14) 


To complete the discussion of the boiling heat transfer we no- 
tice that in the region of single phase gas flow the heat transfer 
is given by equations (3) with appropriately modified param- 
eters. Transition to the single phase heat transfer is initiated 
when dryout transition is detected. 


Further details of the constitutive relations, including pres- 
sure drop correlations, will be provided elsewhere, see 
e.g. (D. G. Luchinsky et al., 2016; D. Luchinsky, Khasin, 
Timucin, Sass, & Brown, 2016). 


3.2. Uncertainties 


Traditionally, the choice of the functional form of the cor- 
relations was the major source of uncertainty in cryogenic 
flow simulations. There have been literally hundreds of 
various modifications proposed for the flow boiling correla- 
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tions (Nellis & Klein, 2009; Shahs, 2006) and the correspond- 
ing functional space is continuously expanding (S. R. Darr 
et al., 2015; Kim & Mudawar, 2014; Konishi & Mudawar, 
2015a). 


In addition, the development of correlations rely on assump- 
tions that can be viewed only as approximations. Further- 
more, these approximations usually do not take into account 
surface wettability, pipe curvature, sub-cooling, and surface 
orientation. 


At the fundamental level the uncertainties in two-fluid mod- 
eling stem from the fact that the interface between two phases 
is continuously fluctuating and neither location nor the shape 
of the interface can be resolved by the model. The inten- 
sity of these fluctuations is especially significant during chill- 
down, when liquid and vapor phases coexist under strongly 
non-equilibrium conditions, see e.g. (Yuan et al., 2007). 


Accordingly, most of the parameters in the correlations dis- 
cussed above should be considered as unknown fitting param- 
eters of the problem. Their relative importance to the sys- 
tem dynamics is also unknown. And a general probabilistic 
framework is required to establish distribution of these pa- 
rameters using available databases of cryogenic flows. 


4. INFERENTIAL FRAMEWORK 


4.1. Probabilistic approach 


The approach to the solution of the inference problem 
adopted in this work is based on probabilistic Bayesian 
method (Ghahramani, 2015). Using this technique one 
can (Ghahramani, 2015) estimate the conditional probability 
of unknown model parameters 6 


P(D|0,m)P(0|m) 


P(6ID,m) = Bm) (15) 
compare different models m 
P(m|D)P 
P(D\m) = ae (16) 


and forecast system response D,, to untested experimental 
conditions 


P(Dy|D,m) = | P(Dyl0,D,m)P(O|D. md. (17) 


In the equations above P(D|6, m) is the likelihood of param- 
eters 0 in model m. In a typical step of parameter inference 
within Bayesian probabilistic framework one has to guess or 
estimate the prior probability of 9 P(@|m) and to obtained the 
improved posterior probability of 9 P(6|D,m) using avail- 
able time-series data D. 


There are two important advantages of this approach. First is 
the method’s ability to select simpler models over more com- 


plex models. Second is the ability to chose model, which is 
flexible enough to resolve the complexity of the experimen- 
tally observed system dynamics. 


4.2. State Space Model 


To perform parameter estimation the equations (1) - (2) are in- 
tegrated over the length of control volumes (D. G. Luchinsky 
et al., 2014c; RELAPS-3D-I, 2012). The result of integration 
at one time step is formally represented in a standard form of 
discrete time state space model (SSM) 


Li41 = f(ti,c) + €t, 
18 
Ue = g(wi, b) +x4- (18) 


Here x; is a set of the main dynamical variables of the sys- 
tem {Pg,pi,Tg, Th, Ug, Ut, D, Ri ag averaged over corre- 
sponding control volume at time instant ¢ and c is the set of 
the model parameters. 


The observations y, in the SSM are related to the unobserved 
states via nonlinear function g(a, b). €, and x4 in equations 
(18) are independent identically distributed sources of Gaus- 
sian noise. 


The simultaneous parameter and state estimations in nonlin- 
ear system (18) can be performed using a number of tech- 
niques, see e.g. (Ghahramani, 2015; Smelyanskiy, Luchin- 
sky, Millonas, & McClintock, 2009; Duggento, Luchinsky, 
Smelyanskiy, & McClintock, 2009). Here, to simplify the 
problem we assume that measurement noise can be neglected 
and that pressure and fluid temperature T; can be measured 
directly in the experiment. 


In the simplest case of general importance the problem was 
reduced to the curve fitting problem, cf (Cullimore, 1998), 
which was solved using global search algorithms available in 


Best: 689365 Mean: 694374 


9 
Number of variables (14) 


Figure 1. (a) Convergence of the genetic algorithm for simul- 
taneous optimization of 14 model parameters. (b) Best values 
of the model parameters. 
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MATLAB optimization toolbox including genetic algorithm, 
pattern search, simulated annealing, and particle swarm. 


An example of convergence is shown in Fig. 1. In this ex- 
ample the genetic algorithm was used for simultaneous op- 
timization over 14 model parameters including characteris- 
tic time of the heat transfer to the wall, scaling coefficients 
for the mass transfer at the interface and at the wetted wall, 
for the film boiling heat transfer, for the dependence of the 
critical heat flux on the void fraction, for the heat transfer 
coefficients on both sides of liquid vapor interface, for the 
heat transfer to the dry wall, for the temperatures of the crit- 
ical heat flux and minimum film boiling, exponent of the 
Reynolds number in the Ditus-Boetler at the dry wall, and 
parameters of the transition boundary to the dispersed flow 
regime. 


To enforce the convergence an extensive preprocessing was 
used including sensitivity analysis of the model with respect 
to the full set of parameters, selection of a subset of the most 
sensitive parameters for each flow regime, and simplified di- 
rect search for initial estimations of sub-optimal values of the 
selected model parameters. Some further details of the ap- 
proach can be found in (D. Luchinsky et al., 2016). 


Below we provide an example of application of this approach 
to the analysis of chilldown in cryogenic horizontal line. 


5. APPLICATION TO THE ANALYSIS OF THE CRYO- 
GENIC LOADING SYSTEM 


The time-series data analyzed in this work were obtained 
at the Simulated Propellant Loading System (SPLS) built 
at NASA KSC (Johnson, Notardonato, Currin, & Orozco- 
Smith, 2012). 


5.1. SPLS 


The SPLS line includes storage (ST) and vehicle (VT) tanks 
connected via cryogenic transfer line that has a number of 
in-line (CV) and bleed (BV) valves that control the flow. The 


T1501 


17202 TT105 cv7 Bv8 
PTI02 PT190 


Figure 2. Sketch of the cryogenic transfer line built at KSC. It 
includes storage tank (ST) and vehicle tank (VT); the in-line 
control valves: CV1 through CV8; remotely controlled bleed 
valves: BV1 through BV8; ten in-line temperature sensors 
(TT) and 9 in-line pressure sensors (PT). 
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Figure 3. Model predictions (green lines) for the fluid temper- 
ature are shown in comparison with the experimental time- 
traces (black lines) at 6 different locations along the transfer 
line. 


flow conditions are monitored using a set of pressure (PT) and 
temperature (TT) sensors. The sketch of the transfer line indi- 
cating the locations of the control valves and sensors is shown 
in the Fig. 2. The typical chilldown sequence measured by 
temperature sensors is illustrated in Fig. 3. The chilldown 
begins at ~800 sec by opening the valve CV1 (note, that the 
initial time is chosen arbitrary). At this time all other valves 
remain closed and after a small temperature drop partial re- 
covery of the temperature measured by sensor TT202 can be 
observed in Fig. 3(a). 


Approximately ~270 sec later the valves CV2 and BD 1, 
2, and 3 are opened allowing liquid nitrogen to flow down 
through the system up to the location of the valve CV3. The 
cooling of the fluid in the pipes can be detected at this stage at 
temperature sensors TT202 though TT774. The chilldown of 
the first part of the transfer line is accomplished ~1000 sec 
from the chilldown initiation, when the temperature sensor 
TT74 indicates the presence of the liquid nitrogen. 


The final stage of chilldown begins at ~2000 sec, when valve 
CV3 is opened and nitrogen is allowed to flow through the 
whole system. The chilldown of the second part of the trans- 
fer line can be observed at the sensors TT 146 through TT226. 
The chilldown is completed when sensor TT226 indicates the 
presence of the liquid nitrogen. 


The comparison of the model predictions (green lines in the 
Fig. 3) with the experimental time traces (black lines) shows 
that the model can quite accurately reproduce all three stages 
of the chilldown. We note that the model integration is fast. 
The integration of nearly 3000 sec of real time shown in the 
figures takes ~10 seconds on laptop. 


A good agreement between model predictions and experi- 
mental time-series data was achieved using a sequence of 
steps as will be discussed in more details below. 
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Figure 4. The simplified direct search algorithm: conver- 
gence of the model predictions (red dashed lines) towards ex- 
perimental data (solid black lines) for gas temperature mea- 
sured at four sensors: (a) TT202, (b) TT105, (c) TT162, and 
(d) TT174. The search is performed for 6 parameters of the 
SPLS model shown in figure (e). The cost function is shown 
in figure (f). (top) Initial state. (bottom) Final state. 


5.2. Optimization of the SPLS parameters 


The initialization and the update of the parameter distribu- 
tions involves the following steps: (i) choice of the model pa- 
rameters; (ii) definition of the objective function; (iii) analysis 
of the initial parameter’s distributions using sensitivity study; 
(iv) estimation of sub-optimal parameter values using simpli- 
fied direct search; (v) global search for the optimal values of 
the model parameters; and (vi) estimation of the variance of 
the optimized parameters. 


There are about 98 parameters in the present model of the 
SPLS. Approximately half of these parameters correspond to 
the correlation relations. The remaining parameters describe 
valve coefficients, minor losses and leaks of the system com- 
ponents. Each parameter is tagged. Any number of parame- 
ters can be selected for optimization by listing corresponding 
tags in arbitrary order. 


The cost function has the form of the sum of square errors on 
discrete time t,, = {to, ..., tw}, cf. (Cullimore, 1998) 


(19) 


where 7; are weighting coefficients, TF »(€) and oe (c) are 
predicted by the model fluid temperature and pressure at k-th 
location and n-th instant of time, pressure and temperature 
with the hat correspond to the measured time-series data, in- 
dex k runs through different locations of the sensors. 


The sensitivity study reveal ~30 sensitive model parameters 
and subsequent optimization was focused on the analysis of 
this subset of parameters. 


The optimization is performed in two steps. At the first step 
we apply simplified direct search algorithm: (i) the parame- 
ter space is discretized on a regular multi-dimensional grid; 
(ii) parameters are optimized one at a time using 1D search; 
(iii) the resulting optimal value of the given parameter is 
used in subsequent searchers; (iv) once the algorithm scanned 
through all the parameters the order of parameters is changed 
and the search is repeated. It was found that this step is essen- 
tial to speed up the convergence and to find global minimum. 
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Figure 5. Comparison of the model predictions for the tem- 
perature time-series data (colored dashed lines) with the ex- 
perimental results (colored solid lines) for the different open- 
ings of the valve cv112: (i) 15% blue line; (ii) 25% black line 
(nominal); (i11) 40% red line. 


The convergence of the algorithm is illustrated in the Fig. 4 
for simultaneous optimization of 6 parameters, which are 
scaling factors for the flow coefficients of the valves: Sm151 - 
input valve CV1; Sc116 for inline valve CV2; Sd120 dump 
valve BV3; Sd117 dump valve BV2; Sd112 dump valve 
BV2; and the minimum film boiling heat transfer - qmfbsc. 


The sub-optimal values of the model parameters found at this 
step are used as seed values for the global search. The global 
optimization was performed using aa number of algorithms 
including genetic algorithm, pattern search, simulated anneal- 
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ing, and particle swarm. An example of convergence of the 
global search for the optimal values of 14 model parameters 
using genetic algorithm is shown in Fig. 1. 


500 1000 1500 2000 0 500 1000 1500 2000 500 1000 1500 2000 


Figure 6. Comparison of the model predictions for the pres- 
sure time-series data (colored dashed lines) with the experi- 
mental results (colored solid lines) for the different openings 
of the valve cv112: (i) 15% blue line; (41) 25% black line 
(nominal); (iii) 40% red line. 


At the final step of the analysis we estimate the variance of the 
optimized model parameters using local search with multiple 
restarts. This completes the analysis of the updated distribu- 
tions of the model parameters given experimental time-series 
data within a general inferential framework. 


This procedure can be systematically continued as soon as 
new experimental data become available. Furthermore, alter- 
native functional forms of the two-phase flow correlations can 
be compared with each other using proposed framework and 
available experimental databases. 


The proposed approach is flexible and allows one to impose 
sufficient number of constraints to guarantee physically valid 
solution of the optimization problem. 


To verify the predictions capabilities of the approach we per- 
formed the following tests. We modify in the experiment 
one of the control parameters while keeping all other controls 
at the fixed values corresponding to the nominal chilldown 
regime considered in the previous section. To make predic- 
tions we keep the model parameters unchanged except the 
one modified in the experiment. 


The results of model predictions are compared with the ex- 
perimental time-series in the Figs. 5 and 6. In this particu- 
lar set of test measurements the opening of the control valve 
CV112 (BV1) was 15% (record 0056) and 40% (record 0057) 
as compared to the nominal opening 25 % (record 0058). 


It can be seen from the figures that the system pressure is 
accurately reproduced for all the time-series data. The model 
predictions for the fluid temperature are less accurate, but are 
still very close to the time-traces observed in the experiment. 


6. CONCLUSIONS 


To summarize, we developed an efficient probabilistic frame- 
work for simultaneous optimization of a large number of cor- 
relation parameters of cryogenic loading system that enables 
accurate prediction of experimental time-series data. The ap- 
proach relies on the fast and robust solver for two-fluid cryo- 
genic flow reported in our earlier work. 


The approach involves the following steps: (i) sensitivity 
analysis of the model parameters, (ii) simplified direct search 
for approximate globally optimal values of these parameters, 
(iii) global stochastic optimization that refines the estimate 
for parameter values obtained at the previous step, and (iv) 
estimation of variance of the model parameters using local 
non-linear optimization. 


We validated this approach by analyzing chilldown in the 
Simulated Propellant Loading System built at NASA KSC. 
We perform simultaneous optimization of a large number of 
the model parameters. We demonstrated the convergence of 
the algorithm towards experimental data at the two main steps 
of optimization and validate its capability to predict “unseen” 
experimental time-traces. 


We note that the proposed approach paves the way to the de- 
velopment of autonomous control and fault management of 
cryogenic two-phase flows in the future space missions. 
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temperature 
pressure 

specific energy 
specific enthalpy 
heat transfer coefficient 
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Reynolds number 
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cross-sectional area 
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Greek 

a gas void fraction 

p density 

T wall shear stress 

T mass flux per unit volume 

Subscript 

g gas 

l liquid 

w wall 

a interface 
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